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Direct search methods for an open problem of 
optimization in systems and control 

Emile Simon and Vincent Wertz 



Abstract 



The motivation of this work is to illustrate the efficiency of some often overlooked alternatives to 
deal with optimization problems in systems and control. In particular, we will consider a problem for 
which an iterative linear matrix inequality algorithm (ILMI) has been proposed recently. As it often 
happens, this algorithm does not have guaranteed global convergence and therefore many methods may 
perform better. We will put forward how some general purpose optimization solvers are more suited 
than the ILMI. This is illustrated with the considered problem and example, but the general observations 
remain valid for many similar situations in the literature. 



Index Terms 

Optimization; Linear systems; Positive filtering 

I. Introduction 

A. Motivation 

There are many open problems of optimization in systems and control theory, often non- 
smooth, non-convex and NP-hard or of unknown complexity. Because of the success obtained 
with linear matrix inequalities (LMIs) since the mid-nineties, a general tendency in the system 
and control literature is to formulate many of these problems with bilinear matrix inequalities 
(BMIs) or LMIs plus a non-convex rank constraint. Once these formulations obtained, the usual 
techniques for trying to solve these problems are iterative linear matrix inequalities algorithms 
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(ILMIs). Such algorithms solve successively different LMI subproblems of the original BMI 
formulation. However, these algorithms loose the guaranteed polynomial time complexity of 
LMI solvers but more importantly most do not have guaranteed convergence to locally optimal 
solutions (i.e. global convergence when starting from any feasible initial solution). 

Such ILMIS, heuristics without global convergence, are then likely to be outperformed by 
other methods like general purpose optimization solvers. This is the fact that we aim at putting 
forward with the current paper: much research work sticks to a unique approach, like BMIs/LMIs, 
and ignores alternative formulations which can however be much more appropriate. 

The performance can be compared in three directions: 1) Objective value, 2) Computational 
time and 3) User time. More details will be given in the paper but we already outline here some 
key elements. 1) The methods that will be proposed are often found in the literature to behave 
well on non- smooth functions and some are proved globally convergent on smooth functions. 
This is not the case for many ILMIs. 2) On non-smooth and non-convex problems there are no 
worst-case complexity bounds neither for the proposed methods nor for ILMIs. Thus, on this 
front, only experience will tell. 3) The time required by the user to implement the method is 
clearly smaller with general purpose optimization solvers, which are intended to minimize any 
user-defined function f{6) : W 1 — > K. Formulating a problem under that form may be done in 
a matter of minutes or hours. But developing a theoretically consistent ILMI algorithm, along 
with a working implementation, can easily require several weeks. This is by the way admissible 
in fundamental research but clearly unacceptable in the industry. 

B. Considered problem 

To illustrate the point of the paper, we will consider a problem for which an ILMI approach 
was recently proposed [12]. This problem is the design of a reduced-order positive filter for 
linear systems. Positive systems, found in many areas (see [0), are dynamic systems with state 
variables and outputs positive at all times. The filter to be designed has to estimate an unmeasured 
output (z) of a system (£) from the measurements (y), with minimum objective value chosen 
as the Hoo norm of the filtering error transfer function (Q) between the exogenous disturbance 
signal (w) and the error (e): the difference between the output (z) and its estimate (z). Details 
on the problem motivations and challenges are given in the introduction of [fT2l . here only the 
elements necessary to reproduce the contribution are recalled. 
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Consider the following asymptotically stable discrete-time system: 



S : < 



x k +i = Ax k + Bw k 
Uk = Cx k + Dw k 
z k = Lx k + Gw k 



(1) 



The notations are classical and identical to those of 02], where all terms are more formally 
detailed. This system is positive if and only if the matrices A, B, C, D, L, G have only positive 
entries flU. Note also that in this paper we only require the positivity of the filter to be designed 
and not of the system S, but it can be assumed because otherwise it is not necessary to design 
a positive filter. 

The aim of the filter is to compute an estimation z\. of the unmeasured signal Zk in S from 
the measured signal y k . More specifically, we want to build the following filter: 

Axk + Byk 

(2) 

Cx k + Dy k 

where A, B, C, D are the filtering parameters to be determined, matrices having only positive 
entries to ensure the positiveness of ik ED, 021. The difficulty of designing this filter stems 
from that particular requirement, which prevents the application of methods normally used for 
this kind of problem (see lfT2l and ref. therein). 



E : 



Xk+l 

z k 



Defining £ A 



xi,x 



TlT 



k ) J 



and ek = Zk — Zk, we get from 0} and © the description of the 



filtering error system: 



where 



A t 



ik+l 

ek 

A 
BC A 
L-DC 



Afik + BfW k 
Cfik + Dfw k 



(3) 



-c 



D 1 



B 

BD 
G-DD 



The transfer function of the filtering error system S e is given by: 

Q{z) = C f (zI-A f )- 1 B f + D f 
The problem considered is defined hereunder. 
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Reduced-Order Minimal Hoc Positive Filtering Problem: 

min 1 1^| |oo s.t. S e stable and A,B,C,D positive (4) 

£ 

This problem is a generalization of the one in lfl2l . The difference is that we do not only seek 
a filter respecting a given level ||£?||oo < 7 but we rather wish to minimize this level. This 
objective is more meaningful and also more practical for the user who does not need to choose 
an arbitrary level 7 beforehand. Note that considering the expression of Af, since S is stable 
we have that the filtering error system S e is stable iff the filter £ is stable. 

II. The general approach 

The problem above can be expressed as either of the two following optimization problems of 
minimizing an objective function f{9) : W l — > R: 
The unconstrained problem: 

{00 if 9 not positive or S e unstable 
\\G\\oo otherwise 

The constrained problem: 

min f{9) = \\G\\oo (00 if S e unstable) s.t. 9 positive 

6 

where n is the number of variables in 9, the vector containing the entries of A, B, C, D. These 
problems can actually be dealt with many general-purpose optimization solvers. Amongst these 
we outline some Direct Search (DS) methods, using the unconstrained formulation: the Mesh 
Adaptive Direct Search (MADS) ED, the Multidirectional Search (MDS) JUL M, the alternative 
directions M, the Nelder-Mead algorithm (NM) [jT4), [0, EH, 023, El, 0, the Particle 
Swarm Optimizer (PSO) [|20l (see ifTOl for the latest survey on DS methods). We also mention 
gradient-based methods, which standard and widely known principles are implemented e.g. 
in the following two functions of Matlab: fminunc for the unconstrained formulation and 
fmincon for the constrained formulation, both methods estimating numerical gradients when 
the expression of the gradient of / is not given (see the Matlab help [fT5l and references therein 
for detailed descriptions). All these methods can be used for the problem discussed here. 

Note that when the expression of an (estimate) of the gradient is known, more evolved 
techniques can be developed (like fl}, 10) to (try to) reduce the computational time, especially 
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for problems with larger number of variables. For complex problems like the one considered, 
the gradient information is useful to accelerate the search but without guarantee to find better 
optima than methods which do not use this information (see benchmarks in IfTTl ). 

A. A method put forward 

For illustration purposes we put forward one method in particular, the Nelder-Mead algorithm 
(NM). This method is probably the most popular algorithm for unconstrained minimization of 
problems of moderate dimension (e.g. < 50 to 100 variables). The central question is whether 
locally optimal solutions can be reached, which is not guaranteed for complex problems with 
most (if not all) iterative Linear Matrix Inequality algorithms (ILMIs) including the one proposed 
in [fT2l . as motivated in the introduction. Note that for further illustration, we will also use several 
other methods than NM in the next section. 

The NM algorithm does not appear often in (recent) systems and control literature. This 
algorithm was first proposed almost fifty years ago in [fl4l and belongs to the class of DS 
methods. These methods, characterized by the fact that they do not build a model of the objective 
function, belong to the broader class of derivative-free optimization methods which do not use 
any gradient or Hessian information. The reader is referred to the last survey on DS methods 
ifTOl Sec. 1.4] for a broader description. The basic ideas behind NM are briefly described in the 
next paragraph, detailed descriptions can be found in [fTTl or 03], lfT6l . 

The first step is the generation of an initial simplex of n + 1 solutions around and including 
the provided initial solution. The objective function is evaluated at each of these n + 1 solutions 
and sorted from the best to the worst. Then NM chooses iteratively between several possibilities 
(or steps) to change the shape of the simplex (eventually displacing it), trying to find better 
solutions. For example the basic step is that the worst solution is reflected on the other side of 
the simplex, in order to create a 'downhill' effect. 

Two easily available implementations of NM are the f minsearch function in the Optimiza- 
tion Toolbox of Matlab 01], 031 or the nmsmax function of [8J. It must be noted however 
that these implementations may fail to converge to locally optimal solutions, starting from a 
feasible initial solution. Indeed the original NM can fail because of the deterioration of the 
simplex geometry or lack of sufficient decrease. However the convergence of the method can be 
guaranteed on smooth functions by taking care of these situations, which is done for example 
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in [fT6l where an additional step ensuring the convergence is proposed. Actually there is a third 
easily available implementation, the grid restrained NM algorithm J3]. This algorithm enjoys, 
like [fT6ll , guaranteed convergence on smooth functions and will also be used in our tests. 

An earlier method to improve the convergence of NM proposed in |[T3l consists of an involved 
restarting strategy. Here we suggest a method with easier implementation: restarting the algorithm 
until the latest obtained objective value is not better than the previous one to a given accuracy. 
This does not formally guarantee the convergence to local solutions on smooth objective functions 
but renders it much more likely. Indeed, restarting NM regenerates its simplex and this allows 
in practice to properly span the space M. n around the last solution. 

Anyway, the hypothesis for convergence analysis in lfl6l . lH3l . [|3l require that the function 
be smooth. And as pointed out in |0Q, the norm objective function is non-smooth and this 
may cause the algorithm -or any of the techniques cited above- to stop at suboptimal solutions. 
This is the reason why we will use not only but also try the idea of restarting NM in the 
results section. 

We summarize that the point is not to guarantee the convergence to locally optimal solutions 
for all optimizations, which require smooth objective functions, but instead to propose an eas- 
ily implementable technique that performs well and in practice often reaches locally optimal 
solutions even for non-smooth objectives functions. 

Nelder-Mead with local restart(s) can be more formally defined as: 

0/(1) =Nelder-Mead(/(#), 6 h options)); 

acc — 1; (> e s ) 

while acc > e s do 
i — i + l; 

6 f (i) =Nelder-Mead(/(0),0/(i - I), options); 

acc = abs(abs(f(6 f (t - l))//(0/(i))) - 1); 
end while 
return 9f(i) 

where 9f(i) is the solution after optimization i, Oi the initial solution, options contains the 
stopping criteria and tolerances of the Nelder-Mead implementation, abs is the absolute value, 
e s is the stopping accuracy required for acc, the relative improvement between the last and the 
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current objective value. 

We already mention another kind of useful restarts, also suggested e.g. in [|T3l . the global 
restarts (or multi-starts). These consist simply of running other optimizations from other different 
initial solutions 0j. The aim of these global restarts allows to better span the search space which 
is useful to try to: 1) get around the non- smoothnesses that could have caused other optimizations 
to stall, for non-smooth objective functions 2) to find better local optima (nearby), for multimodal 
objective functions 3) get away from suboptimal solutions against the infinite penalization barrier, 
which is an implicit constraint. By obtaining an increasing number of final solutions having the 
best objective value found so far (to a given accuracy), we get increasing probabilities that these 
solutions are locally optimal and in a lesser measure globally optimal. 

B. The advantages 

This illustrating algorithm may compete with other methods, e.g. with [6] using benchmarks 
of static output feedback optimization in ifTTll . It can often find the same objective values or even 
better ones, possibly in shorter computational times although this is more an exception rather 
than the rule, the main drawback being the lack of explicit handling of non- smoothnesses as 
proposed in [HI, [0. The main advantage is the great flexibility to handle any objective function 
f(0), only some of which can be written under other frameworks like LMIs. Moreover, even 
without using gradient information, it can often lead to locally optimal solutions. Also, unlike 
usual quasi-Newton methods, NM has the ability to explore neighboring valleys with better local 
optima and likewise this exploring feature may allow NM to overcome non-smoothnesses. 

The flexibility is what makes it a candidate of choice for dealing with designs requiring 
particular structures or properties, such as ensuring the positivity of the solution like desired 
in this paper, while exploring a non-smooth objective function. More precisely, NM will be 
particularly efficient for optimization objectives as an alternative to ILMIs, as motivated before 
and described further as follows. 

ILMIs are often efficient for feasibility problems such as finding a stabilizing controller, 
possibly also ensuring a given performance level, even though they may fail to find a solution 
even if there exists one (see flU, flU). However when used for minimization problems, like 
minimizing a norm of a performance channel, ILMIs have in general no guarantee of convergence 
to locally optimal solutions (as discussed in [fT8ll and noted in [Q21). For instance in [12] this 
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stems essentially from the fact that the objective minimization would require a combinatorial 
optimization problem involving the objective (7) and a parameter (e). Considering this important 
drawback, other methods should be used instead of ILMIs for objective minimizations. ILMIs 
can be useful when obtaining feasible (e.g. stable) initial solutions is not trivial, which is not 
the case here (the system is already stable). And then more suited methods should be used for 
optimizing the objective, like DS methods. 

Indeed using these methods can often lead to locally optimal solutions, unlike the algorithm 
of lfl2ll that produces suboptimal solutions under an a-priori defined level 7. Although the non- 
smoothness prevents guarantees of local optimality for direct search methods, it is conjectured at 
places in the literature that this does not happen often (see ref. in 0]). How 'often' depends on 
the density of non-smoothnesses in the objective function and in particular at 'partial' optimal 
solutions, optimal in some but not all directions. This is actually a matter of debate, on which the 
reader is advised to consult e.g. JT), [0 and the references therein. And as noted before, global 
restarts can be used to try to reduce the likeliness of getting stuck at a suboptimal solution. 

Many optimization problems in systems and control can be quickly implemented to be solved 
with DS methods. Indeed the objective function is easily built and evaluated with adequate 
methods, for example the functions from the Control System Toolbox of Matlab. So there is no 
need to modify or add cumbersome LMIs, which are also especially large with the technique of 
|fl2l (see the matrix inequality (16) there). Likewise the unhandy implementation of an ILMI is 
avoided, moreover often involving parameter(s) for which no automatic a priori choice can be 
made. In short, the method can be implemented in a brief time even by non-expert users. Also 
DS methods are useful for simulation-based optimization [[TOl Sec. 1.21], where more intricate 
objective functions that could hardly be obtained from analytical expressions are computed using 
simulation results (illustration given in next section). 

An important remark can be made regarding the worst-case computational complexity bounds. 
In [H we can read that MDS and NM are not competitive with more sophisticated methods such 
as (quasi-)Newton methods when applied to smooth problems. Indeed it seems clear that gradient 
expressions and moreover Hessian information should accelerate the search. However a recent 
result [ETI proves that directional DS methods share the worst-case complexity bound of steepest 
descent for the unconstrained minimization of a smooth function. This gives further motivation 
to consider DS methods. 
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Let us also mention that the formulation of the problem in [[121 requires many matrix variables 
Pi, P2, R, U, V, Pj, i — 1, 6 and three scalars a, (3, 7 -with U, V, a, (3 entering non-affmely 
the matrix inequalities, thus the need of an ILMI- whereas here we only need the original 
variables A,B,C,D which put together lead to the same size as only R in |[T2l . Using LMIs 
leads to this typical key problem of inflation of size and number of variables for large systems 
(E and E). Instead, with general purpose solvers, only the original problem's variables (E) are 
used. 

LMI problems are solved efficiently with interior point methods converging to the optimum 
in worst-case polynomial time. This is not guaranteed with the method proposed here, but also 
not with ILMIs, since the iterative scheme does not preserve the polynomial time complexity. In 
practice NM deals well with problems with limited number of variables (e.g. < 10-20) but its 
performances decrease notably for larger number of variables (e.g. > 50-100). As noted before, 
one should then use the advantage that DS methods can handle problems formulations with the 
least number of variables. 



III. Results on example 
The positive system considered is given by 02): 



Xk+l 



0.1595 0.1890 0.2713 



0.5091 











0.6740 







1 




Vk = 




Xk + 




1 





Xk + 



0.1350 0.0128 
0.3850 1 0.0510 
0.1021 0.1250 



0.1250 
0.1460 



The objective is to estimate z k = [1 0]x k = x\. For the problem considered, any positive and 
stable filter can be used as initial solution. Since like in [12] we want to design a first order filter, 
we can choose the scalar A in [0, 1[ and the other variables in M + . We then simply use rand 
of Matlab to generate the random initial solutions (entries chosen uniformly in [2~ 53 , 1 — 2~ 53 ]). 
Performing three optimizations, each from a different random initial solution, gives us the best 
following solution: 



'This value is erroneously written 0.0128 in 11121 . 
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A = 0.06978, B = [0.53667 2.13004] , C = 0.15218, D = [0.15435 0.10931] 

Instead of the solution in Ifl2ll : 

A = 0.22819, B = [0.00003 0.00003] , C = 0.14130, D = [0.17889 0.34404] 

Let it be clear that we do not need the solution of lfi~2l as initial solution, indeed here we simply 
use independent random feasible solutions and therefore the proposed technique is self-standing. 
Perhaps for some particular problems ILMIs would be necessary to find initial solutions, but 
this is not the case here. We nevertheless give as further information what happens when using 
the solution of |fi~2| as initial solution, at the end of this section. 

Note also that since linear systems admit an infinity of equivalent state-space representations, 
we can restrict the search space to one representation (e.g. by fixing here C = 1). On the one 
hand, reducing the number of variables reduces the computational time. On the other hand, this 
removes the possibility for DS methods to explore different state-space coordinates, therefore 
reducing their ability to find better local optima and getting around non- smoothnesses or away 
from the border of (implicit) constraints. 

Two of the three solutions found have a Hoc performance level around 0.0447, significantly 
better than the level 0.1415 reached in [fi~2l . We present in Fig. [TJ the same simulation as in 
021 of the actual state x\ and its estimations, using the following initial condition of the error 
system: [0.03, 0.08, 0.10, 0.05] T and exogenous disturbance input: w k = [1/(1+0.25*;), e -° 02fc ] T - 

The gain of performance level can be seen, with a maximum absolute error about three 
times smaller with the restarted NM. A third simulation response is also provided, obtained 
by minimizing an objective specific to this simulation: the sum of the absolute values of the 
errors e k . The obtained filter is given by: A = 0.011, B = [0.00598 1.779], C = 0.3, D = 
[0.000046 0.0607]. This gives an illustration of simulation-based optimization. 

To give further illustration of the performance of the proposed method, we run 100 optimiza- 
tions each starting from a different random initial solution. We then give the minimum and the 
average of all 100 obtained objectives values. Also we give the percentage of objectives values 
that were smaller than 0.1415 and of those smaller than 0.0448. As last indication, we give the 



January 20, 2013 



DRAFT 



11 



0.5 1 1 1 1 i i i i i i 




Time k 

Fig. 1. State x k , estimations and corresponding errors e^. 



average of the 100 computational times required in seconds^. We perform these tests with the 
proposed technique as well as with other optimizations techniques, each using the same 100 
random initial points. The results are given in Table [I] 



TABLE I 

Comparison of objective values and computational times 



Method 


Obj.: min/mean 


<0.1415 


<0.0448 


Time (s) 


1) 


fminsearch 


0.04470897/0.07 


96% 


5% 


9.57 


2) 


fminsearch + restart(s) 


0.04470746/0.04559 


100% 


73% 


44.2 


3) 


nmsmax 


0.04470746/0.0514 


98% 


38% 


8.53 


4) 


nmsmax + restart(s) 


0.04470746/0.04471 


100% 


100% 


41.7 


5) 


mdsmax 


0.04471485/0.0479 


100% 


5% 


6.81 


6) 


gridnm 


0.04470746/0.045 


100% 


79% 


22.8 


V) 


fminunc 


0.04473497/0.272 


56% 


1% 


2.41 


8) 


fmincon 


0.04470795/0.0627 


98% 


38% 


2.37 



The methods are named after their Matlab implementation file: 1) fminsearch, the Matlab 

'The computer used is a HP Compaq dc7800©, processor Intel Q9300©, 2.5GHz, 3.48Go RAM, software MATLAB 2007b© 
January 20, 2013 DRAFT 
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Optimization Toolbox NM implementation 2) The same with local restart(s), as described in 
the previous section 3) nmsmax, N. Higham NM implementation [8] 4) The same with local 
restart(s) 5) mdsmax N. Higham MDS implementation |[8]| Ifl9ll 6) The grid restrained NM 
algorithm gridnm [3|, proved convergent on smooth functions 7) fminunc also using the 
unconstrained formulation (using 10 100 instead of oo for constraint penalization) 8) fmincon, 
the only method using the constrained formulation of the problem. 

All accuracies required: 'TolF', 'TolX', 'TolCon' (see Matlab help [ T51 for descriptions), e s and 
the accuracies in nmsmax, mdsmax, gridnm have been set to 1CT 7 as well as the tolerance 
of the objective ||{?||oo (norm(G,inf,tol)) evaluation. Also the numbers of iterations and 
functions evaluations were not limited (e.g. 'Maxlter' and 'MaxFunEvals' set to inf). 

Except for 7), at least 96% of the trials give a better performance level than 0.1415 from 
021. This is a first illustration of how general purpose optimization solvers might compete with 
ILMIs. Observe then the results 2), 4) and 6), obtained with improved NM. As can be seen 
these have high success rate at reaching objectives values very close to the best one found. 4) 
is especially good since it was 100 % successful. Considering this, even though f(9) is non- 
smooth here, the solutions found with objective value < 0.0448 are very likely locally optimal 
and probably globally as well. This illustrates the efficiency of the method put forward. 

The average computational times are under the minute, therefore the methods are not only 
convenient to encode but also reasonably fast to run. Note that the computational times could be 
shortened by using gradient expressions like in [fl], [ID or, as suggested in HI, by interrupting the 
bisection algorithm computing the norm once it is sure that the solution being evaluated is 
better or worse than the other solutions. According to one of the authors in [fl2l the computational 
time required there was around a few seconds, which is comparable to those required here. 

We also run the different alternatives starting from the solution in |[T2ll . which gives the 
following improved objectives values: 1) 0.1395 2) 0.0565 3) 0.0455 4) 0.0447 5) 0.0466 6) 
0.0449 7) 0.1391 8) 0.0454. We see that two methods were unable to improve the solution of 
|[T2l to be close to the objective value 0.0447, which gives an indication that this solution is 
badly located. Indeed its very small B matrix, 'blocked' near the positivity constraint, almost 
cancels the effect of the dynamical part of the filter: the solution in [12J can be approximated by 
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a static filter (the D matrix) with almost no impact (0.003%) on the objectivaj. On this matter 
we also mention that an important element for the performance of ILMIs is to try to stay away 
from the borders of the BMI feasible set (using e.g. |[T8l Sec. 4.3]), otherwise the algorithm may 
stall early at partial optimal solutions. 

IV. Conclusion 

Considering the heuristical nature of ILMIs, which can be efficient for feasibility problems 
but have in general no guarantee of convergence to locally optimal solution for minimization 
objectives nor a bound on the worst-case complexity, many other methods can be used as 
competitors of such techniques. 

The problem of [12], designing a reduced-order positive filter to estimate the output of a 
positive system under a given maximum error level 7, is dealt with in that paper using an 
ILMI and indeed it can be read there that the convergence of the algorithm is not guaranteed. This 
motivates us to propose techniques where the error level is minimized. We then consider 
a less restrictive approach than (I)LMIs and reformulate the problem so as to be solvable by 
general purpose optimization solvers. In particular, we put forward the NM algorithm improved 
with local restarts that can behave well even with non-smooth objective functions. 

This approach has apparently only advantages compared to that of [12]. It is easy to encode 
and use even by non-expert users. It is more flexible, i.e. straightforward to modify for example 
to 1) handle continuous-time systems 2) change the objective into H 2 minimization or multi- 
objectives 3) take into account complex requirements such as a structure of the solution. The 
inflation of size of the system -but not necessarily of the filter- will have a smaller impact 
on the computational time than techniques using LMIs, since the only variables are those of 
the filter to be designed (no additional variables needed, like Lyapunov matrices). Also the 
technique often leads to locally optimal solutions, where 'often' depends on the objective to be 
minimized, whereas this should be seldom the case with most ILMIs (see |[T8l and references 
therein, in particular [7]). Finally since it is generally fast for moderate size problems, it can be 
run multiple times from several feasible initial solutions. And so by getting several times the 
same best solution one gets an increasing probability that this solution is indeed locally optimal 
or even globally optimal. 

2 With value actually around 0.1417 and not 0.1415, using the tolerance 10 -7 in norm instead of the default 10 -2 
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In the end our claim is that most open optimization problems in systems and control should 
not be handled with conservative LMI formulations, except for providing an initial suboptimal 
solution. Indeed conservative LMI formulations or ILMIs to deal with BMI problems typically 
do not have guaranteed convergence to local optima. To get locally optimal solutions one will 
need to use other approaches, the best two methods being probably 0], 10, which use gradient 
expressions and take non-smoothnesses into account. Nevertheless the method proposed can be a 
competitor to these techniques by finding the same or better objectives. Also the set of problems 
that can be solved by DS methods is broader than those dealt with the methods in JT), [0, thanks 
to the great flexibility offered by the general objective formulation f{6) and without need of 
gradient expression. 
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